4j 


JIM AM 


ZlZlKp 


1028724 


AIAA-2002-0383 

LEWICE 2.2 Capabilities and Thermal Validation 

William B. Wright 
QSS Group, Inc. 

Cleveland, OH 


40th Aerospace Sciences 
Meeting & Exhibit 

January 8-1 1 , 2002 / Reno, NV 


For permission to copy or republish, contact the American Institute of Aeronautics and Astronautics 
1801 Alexander Bell Drive, Suite 500, Reston, VA 20191-4344 


AIAA-2002-0383 


LEWICE 2.2 Capabilities and Thermal Validation 

William B. Wright* 

QSS Group, Inc. 


Cleveland, 

I. Abstract 

A computational model of bleed air anti-icing and 
electrothermal de-icing have been added to the 

LEWICE 2.0 * 1 software by integrating the capabilities 
of two previous programs, ANTICE 2 and LEWICE/ 

Thermal 3 . This combined model has been released 
as LEWICE version 2.2. Several advancements have 
also been added to the previous capabilities of each 
module. This report will present the capabilities of the 
software package and provide results for both bleed 
air and electrothermal cases. A comprehensive vali- 
dation effort has also been performed to compare the 
predictions to an existing electrothermal database. A 
quantitative comparison shows that for deicing cases, 
the average difference is 9.4°F (26%) compared to 
3°F for the experimental data while for evaporative 
cases the average difference is 27°F (32%) com- 
pared to an experimental error of 4°F. 

II. Nomenclature 

A Area(m 2 ) 

C e Mass concentration of water in air at the 

edge of the boundary layer (kg/m 3 ) 

C p Heat capacity (kJ/kgK) 

C s Mass concentration of water in air at the 
surface (kg/m 3 ) 

e Evaporative pressure (N/m 2 ) 

F Force per unit span (N/m) 

H Enthalpy per unit volume (kJ/m 3 ) 

h Heat transfer coefficient (kW/m 2o K) 

h m Mass transfer coefficient (m/s) 

k Thermal conductivity (kW/m°K) 

L Lewis number 

Lf Latent heat of fusion (kJ/kg) 

L v Latent heat of vaporization (kJ/kg) 

LWC Liquid water content of air (kg/m 3 ) 

M Mach number 

rit Mass flux (kg/m 2 s) 


Ohio 44135 


MW Molecular weight (kg/kg mol) 

N Mass flux (kg/m 2 s) 

N f Freezing fraction 

q” Heat flux (kW/m 2 ) 

q’” Volumetric heat source (kW/m 3 ) 

Pr Prandtl number 

R Ideal gas constant (kJ/kg mol°K) 

r Recovery factor 

T Temperature (°K) 

t Time (s) 

V Velocity (m/s) 

W e Weber number 

x Parallel to the airfoil surface (m) 
y Perpendicular to the airfoil surface (m) 

Greek Letters 

Y Ratio of heat capacities Cp/C v 

p Density (kg/m 3 ) 

Subscripts 

a accretion 

aero aerodynamic 

air air 

conv convection 

e edge of boundary layer 

evap evaporation 

fh frictional heat 

ice ice 

k layer number 

lat latent heat 

I liquid phase 

m melt phase 

nc net amount of convection 

o total property 

r melt range 

rec recovery property 

s solid phase 

w water 

x x-dependent 

y y-dependent 

°o free-stream property 
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III. Introduction 

The removal and/or prevention of ice on aircraft 
components is vital to aircraft performance and oper- 
ation. Even small amounts of ice can have disastrous 
consequences. Because of this, several methods of 
ice prevention and removal have been designed. 
Methods of ice control can be arranged into two 
broad categories: anti-icing methods and de-icing 
methods. Anti-icing methods are concerned with the 
prevention or minimization of ice buildup on pro- 
tected surfaces. De-icing methods are concerned 
with ice removal after and during ice build up. 

The primary means of preventing ice formation 
on wings and engine inlets for modern commercial 
transport aircraft is by extracting hot air from the 
compressor and blowing it on the inside surface of 
the leading edge through small holes drilled in a pipe. 
The amount of heat supplied is usually determined 
by correlations based upon hole and pipe diameter, 
the number of holes, the mass flow rate of the air and 
the supply temperature. Ice accretion is prevented by 
supplying enough energy to evaporate the impinging 
water (evaporative anti-icing) or by maintaining a sur- 
face temperature above freezing (running wet anti- 
icing). 

Another widely used method for either anti-icing 
or de-icing aircraft components is with an electrother- 
mal pad. By this method, heater mats are installed 
beneath the skin of a wing surface surrounding the 
leading edge as shown in Fig. 1 . Thermal energy in 
the form of conducted heat melts a thin layer of ice 
which destroys the adhesion force at the ice-surface 
interface. Aerodynamic forces then sweep the ice 
from the surface. When more heat is supplied, the ice 
will melt completely and runback to unheated 
regions. If enough heat is supplied, the water will not 
freeze on the surface, creating an anti-icing condition 
similar to the bleed air case. 

The report is divided into three sections. The first 
section will provide a description of LEWICE 2.2 
model, with emphasis on the anti-icing and de-icing 
physics. The second section will provide several 
examples illustrating the capabilities of the bleed air 
anti-icing model, including design analysis performed 
for an IRT test scheduled in FY02. The third section 
will provide example cases and comparisons from 
the electrothermal model. A detailed comparison with 
experimental data is performed, an effort comparable 


in complexity and scale to the validation effort per- 
formed for LEWICE 2.0 4 . 

IV. LEWICE 2.2 

The computer program LEWICE embodies an 
analytical ice accretion model that evaluates the ther- 
modynamics of the freezing process that occurs 
when supercooled droplets impinge on a body. The 
atmospheric parameters of temperature, pressure, 
and velocity, and the meteorological parameters of 
liquid water content (LWC), droplet diameter, and rel- 
ative humidity are specified and used to determine 
the shape of the ice accretion. The surface of the 
clean (un-iced) geometry is defined by segments 
joining a set of discrete body coordinates. The soft- 
ware consists of four major modules. They are 1) the 
flow field calculation, 2) the particle trajectory and 
impingement calculation, 3) the thermodynamic and 
ice growth calculation, and 4) the modification of the 
current geometry by addition of the ice growth. 

LEWICE applies a time-stepping procedure to 
“grow” the ice accretion. Initially, the flow field and 
droplet impingement characteristics are determined 
for the clean geometry. The ice growth rate on each 
segment defining the surface is then determined by 
applying the thermodynamic model. When a time 
increment is specified, this growth rate can be inter- 
preted as an ice thickness and the body coordinates 
are adjusted to account for the accreted ice. This pro- 
cedure is repeated, beginning with the calculation of 
the flow field about the iced geometry, then continued 
until the desired icing time has been reached. 

The thermal module calculates the 2D transient 
(time-dependant) heat transfer in a body. It can han- 
dle multiple composite layers, where each layer can 
have different thermal properties including tempera- 
ture dependence and anisotropy. Each section can 
contain a heater which is individually controlled with 
separate on/off times and power densities. Ice 
growth can be predicted with or without heaters, and 
the model includes various modes for ice shedding 
and water runback. It can also function as an ice 
accretion simulator with heater power turned off. The 
original deicer program has been documented in pre- 
vious reports 3,5 . For version 2.2, this module has 
been modified to include boundary conditions which 
model the heat flux or heat transfer coefficient from a 
bleed air system. 
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LEWICE 2.2 has been designed to combine the 
best features of LEWICE 2.0, LEWICEAThermal and 
ANTICE. It incorporates the multi-time stepping 
capabilities of LEWICE as well as the ability to easily 
interface with other software. It contains the flexibility 
in heater design and materials from LEWICE/Ther- 
mal and incorporates the bleed air analysis and rivu- 
let runback capabilities found in ANTICE. 


foil with electrothermal heating can be represented 
as: 


, _ . 3T , 3T . d T m 
p p )k at x k 3x 2 +k y’ k a y 2 +q k 


(i) 


For the ice layer, the governing equation in terms 
of the enthalpy is given by 


V. Thermal model description 

The following assumptions were made in the 
development of a mathematical model for heat con- 
duction in a composite airfoil: 

1 . The thermal physical properties of the mate- 
rial composing each layer of the airfoil may be differ- 
ent. 

2. Each layer is either isotropic or orthotropic, 
meaning that cross-derivative terms in the anisotro- 
pic heat equation are neglected, but that thermal 
conductivity with the grain can be different from the 
thermal conductivity against the grain of the material. 

3. There is perfect thermal contact between lay- 
ers. 

4. The thickness of each layer can be assumed 
constant in the chordwise direction. 

5. Curvature effects of the airfoil are not taken 
into account because the deicer thickness normal to 
the airfoil is thin compared to the effective airfoil thick- 
ness. 

This assumption does not mean that curvature 
effects are unimportant nor that curvature effects are 
ignored. The assumption is that body curvature only 
effects the air flow, the collection efficiency, the con- 
vective heat transfer and the ice shape. The conduc- 
tion heat transfer can be modeled as a two- 
dimensional rectangular entity if the deicer is thin 
compared to the airfoil. This is usually true, since a 
thick deicer is not very efficient. 

6. Heat transfer in the spanwise direction is neg- 
ligible. 

7. The ambient conditions (temperature, velocity, 
etc) are constant with respect to time. 

8. The phase change of the ice is assumed to 
occur over a small temperature interval near the true 
melting point rather than at the melting point itself. 

With the above assumptions, the mathematical 
formulation for the problem of unsteady heat conduc- 
tion in a chord-wise two-dimensional composite air- 


3h , faV ajr) 

at ice {dx 2 + a y 2 J 


( 2 ) 


In order to solve this problem, it is necessary to 
remove enthalpy from the governing equation above 
and replace it with temperature. The standard rela- 
tionship between enthalpy and temperature is 



(pCp)J 

T<T m 

H = 

(pC p ) s T m < H < PiC ps T m + p,L f 

T m 


_(PCp) s T m + PiL f + (pC p ),(T-T m ) 

J>T m 


However, this relationship is nonlinear, as there 
are multiple values for the enthalpy while the temper- 
ature remains at the melting temperature. The 
numerical method employed requires that a linear 
relationship exists, so that the coefficient matrix cre- 
ated can be easily inverted. To accomplish this, ice is 
assumed to melt over a very small temperature range 
near the melting temperature instead of at the melt- 
ing temperature. The modified form of the enthalpy/ 
temperature relationship is: 


(pC„> t 

P S 

< T -T m > 

p m C p s T m + p w L f T f 
p / C p s T m + p l L f + ( pC p\ (T - T m -V 


T<T 

m 

T < T < T + T 
m m i 

T > T + T 
m r 


(4) 


It is convenient to define a specific heat capacity 
in this melt region, which is given by: 


(c P ) m 


(Pi~Ps)Cp s T m + P|L f 

PmT r 


(5) 


where p m is the density of the region between 
solid and liquid, henceforth to be referred to as the 
‘melt’ region. Therefore, the governing equation for 
the ice layer can be written as: 
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£<p. c n., T > 


-It 

dx‘ ce dx dy lce 



where Eq. 8 or Eq. 9 is used depending upon 
whether the interface in question is in the x-direction 
or in the y-direction. 

Fixed temperature interfaces are given by: 

^surface = T - ( 10 > 

(6) or 


0 

[(pCp) m -(pC p ) b ]r m 

|(pc p ),- ( P c p ) s ]r„, + [( P c p)/ -(pc p )J r, 

where the top expression in each bracket is used 
for the solid phase, the middle expression is used for 
the melt phase, and the bottom expression is used 
for the liquid phase. As long as this range is small, 
the accuracy of the solution is not significantly 

altered. A melting range of 10' 4 K or smaller has 
proven to be necessary for this application. Since the 
temperature is not known prior to this calculation, the 
phase of the ice (solid, melt, or liquid) must first be 
assumed, the temperature calculated using the 
appropriate equation, and the phase checked at the 
end. This creates an iterative scheme to solve for the 
temperature within the composite body. 


Boundary Conditions 

There are four types of boundary conditions 
used for the composite body. They are: 

1 . Interior interfaces between layers of the com- 
posite body. 

2. Fixed temperature outer boundary conditions. 

3. Convection outer boundary conditions. 

4. Heat flux outer boundary conditions. 

At interior interfaces, the temperatures and heat 
fluxes are continuous, i.e.: 


T) 


layerl 


-tin 


-k 


X dx J 

dr 
y dv . 


layer \ 


layer 1 


^)layer2 

(7) 

k ^ T ] 


'dx ) layer! 

(8) 

dr\ 


'dy /layer! 

0) 


T) surface = ? rec <"> 

Equation 10 is rarely used for the simulation of a 
de-icer pad, but has been incorporated into the pro- 
gram which has been developed so that other heat 
transfer phenomena may be studied by the user. 
Equation 1 1 can be used to model a dry air boundary 
condition. 

Convection boundary conditions are given by: 

-fcjn = h(T,-T m ) (12) 

OX / surface 

= h(T s -TJ ( 13 ) 

OV J surface 

Insulated boundary conditions are obtained by 
setting the heat transfer coefficient in the above 
equations to zero. Additionally, the heat transfer coef- 
ficient can be constant or more likely a function of 
position on both the external boundary and on the 
internal boundary when analyzing a bleed air case. 

Heat flux boundary conditions are given by: 


, an 

X dx ) surface ^ 

an 

y dy) surface 


The heat flux can also be given as a function of 
position as would be the case when performing a 
bleed air analysis. The final type of boundary condi- 
tion used at the top surface of the ice, and it is 
derived in the next section. 

Accretion Boundary Condition 

In early numerical models of a de-icer pad, the 
outer boundary condition at the top surface of the ice 
was assumed to be given by Eq. 12 of the previous 
section. However, this is not adequate when the 
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actual physics of the flow at the exposed surface are 
considered. The approach used here is a modified 
version of the approach used in the LEWICE model 1 
and by Gent and Cansdale 6 . There are several 
assumptions which need to be made in the develop- 
ment of an improved icing boundary condition. These 
are listed below. 

1. The terms in the energy equation which are 
considered important are: 

a. convection losses. 

b. kinetic heating. 

c. evaporative/sublimative cooling. 

d. sensible heat gain/loss. 

e. viscous losses. 

f. latent heat gain. 

g. conduction into or out of the body 

All other terms are neglected. 

2. Kinetic heating is adiabatic. 

3. The Chilton-Colburn analogy is used to relate 
the mass loss caused by evaporation to a heat loss 
caused by evaporation. 

4. Viscous dissipation is given by a ‘recovery fac- 
tor’ as defined by Schlichting 7 . 

5. Air is a perfect gas. 

6. The Antoine equation is used to relate vapor 
pressure to static temperature. This implies that the 
molar volume of the liquid is negligible compared to 
that of the gas and that the heat of evaporation is 
independent of temperature. 

7. Physical properties of air and water, except for 
air density, are independent of temperature. 

8. Mass transfer is proportional to the concentra- 
tion difference across the boundary layer. 

9. The ambient conditions of velocity and mass 
density are constants. 

If dry air flows over the airfoil, the surface tem- 
perature would be increased beyond that of the ambi- 
ent because of the compressible flow of air which 
strikes the surface of the airfoil, thereby imparting the 
kinetic energy of the air to the airfoil. Since this is pre- 
sumed to be an adiabatic process, the local static 
temperature at the edge of the boundary layer would 
be given by: 


T. = 




(16) 


Schlichting relates this temperature to the sur- 
face temperature through the ‘recovery factor’ r which 

is equal to Pr 1/3 for a turbulent boundary layer and 
Pr 1/2 for a laminar boundary layer. This is presumed 
to account for the viscous dissipation in the boundary 
layer. The final result is called the recovery tempera- 
ture and is given by: 


r.(i*rV«;) 

(I+V"') 


(17) 


where T rec is the recovery temperature. The 
amount of energy increase owing to the combined 
effects of the frictional heating by air and the viscous 
dissipation of the boundary layer is normally placed 
in terms of a convective flux: 


q" fh = h(T rec -TJ (18) 

The convective flux is described by Newton’s 
Law of Cooling: 


q" conv = h{T s -TJ (19) 

When convection is combined with frictional 

heating, Eqs. 18 and 19 may be united into a net con- 
vective flux which is given by: 

q" nc = h(T s -T rec ) (20) 

An additional kinetic heating term is also present 

owing to the impact of the water droplets on the sur- 
face. This is given by: 


The heat loss caused by evaporation is deter- 
mined by finding the mass transfer rate of water 
vapor leaving the surface of the airfoil and multiplying 
it by the heat of vaporization. Sublimation is handled 
in the same manner, but using the heat of sublima- 
tion. Since the latent heats of vaporization and subli- 
mation are large, it takes only a small amount of 
mass to be removed by either process before a sig- 
nificant drop in surface temperature is realized. The 
driving force for this process is the concentration dif- 
ference across the boundary layer, C s -C e , where the 
mass density, C, is given by: 
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c = 


e(MW) w 

RT 


( 22 ) 


The evaporative mass flux of water vapor 
through the boundary layer can be described by: 


ing of the supercooled liquid up to freezing, the freez- 
ing of the water into ice (latent heat), and the 
subsequent cooling to the surface temperature. 
These three terms can be expressed in equation 
form as: 


N = h m (C s -C e ) (23) 

The heat loss is given by: 

<7 "evap = L V N < 24 > 

Combining Eqs. 22 through 24, the evaporative 
heat loss can be written as: 


= L v h n ,(MWU e s 

<7 evap ft ^ J J J 


(25) 


The mass transfer coefficient can be replaced 
with the heat transfer coefficient via the Chilton-Col- 
burn analogy 8 between heat and mass transfer: 


h„, = 


(pCp) a j r L 


2/3 


(26) 


Replacing the density of air via the ideal gas law 
gives the net heat lost by evaporation as: 


q"sens I = pl (T m -TJ (29) 

<l"la, = (30) 

q'sensl = ->»u C ps(T s ~ T m ) ( 31 ) 

If the surface is at the freezing point and hence 
only part of the ice freezes, the terms are the heating 
of the supercooled liquid up to freezing and the par- 
tial freezing of the water into ice (latent heat). These 
terms can be described by: 

q”sens3 = « C pl( T n, ~ T J (32) 

<7 "la, = N f< L f ( 33 ) 

Since the assumption has been made that ice 
melts over a small temperature range, the freezing 
fraction can be written in terms of the surface tem- 
perature: 


L v h(MW) w 


evap 


[P 0 (MW)C p ] ai L 


273(^7- e e) 


(27) 


The evaporative pressure at any location can be 
evaluated solely in terms of the local static tempera- 
ture using Antoine’s equation: 


log € s = 4 - JL. (28) 

where A, B, and C are empirical constants. Since 
the temperature at the edge of the boundary layer is 
given by Eq. 13 and the evaporative pressure is given 
by Eq. 25, the only unknown in Eq. 27 is the surface 
temperature. 

The sensible heat transfer and the latent heat 
transfer are determined by tracing the thermody- 
namic path the incoming liquid takes to get to the sur- 
face temperature. The exact form of these equations 
depends upon how much heat is available after the 
other heat losses/gains have been considered. If the 
surface temperature is below the freezing point, the 
sensible and latent heat terms needed are the heat- 


T„, + T - T. 

N _ in I v 

f T. 


(34) 


Hence, the latent heat term for this case can be 
written as: 


( T„. + T T,)m w L f 

„ _ v m i ii i 

*7 lai — -r ( J5 J 

1 r 

If the surface temperature is above freezing, the 
sensible heat is given by: 

4'sensl = -m w C pl (T s -T m ) (36) 

Additionally, there is a sensible heat transfer 
between the runback water and the region into which 
the water flows. This heat exchange is caused by the 
temperature distribution along the airfoil. These 
terms cannot be discounted, as they are quite signifi- 
cant for cases involving a considerable amount of 
runback. The form of the equations describing this 
sensible heat transfer is the same as those for the 
sensible heat terms above, but replacing m w with the 

mass flux of runback water, m rb , and the ambient 
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temperature, T Ml with the temperature of the runback 
water, T^. 

The accretion boundary condition can be found 
by setting the heat flux at the outer surface equal to 
the sum of the heat fluxes described above. This can 
be written as: 

_ (^' 9 5 ) — ^ nc evap ~ q ke~ Q lat — tf sens 1 ^) 

Note that the heat flux is directed outward, such 
that a term which results in a heat gain at the surface 
has a negative sign and a term which results in a 
heat loss has a positive sign. As with the governing 
equations for conduction within the de-icer, the sur- 
face temperature is not known prior to this calcula- 
tion, so that the Method of Assumed States must be 
used here as well. 

Accretion Mass Balance 

In the previous section, conservation of energy 
was applied at the top surface of the ice. This 
resulted in the derivation of an appropriate boundary 
condition for the heat conduction equation. In this 
section, conservation of mass will be applied to the 
same surface to determine the amount of ice which 
forms on the airfoil. The additional assumptions 
which have been made in this derivation are: 

1. The flow of surface water is shear driven. 
Hence if all of the incoming water does not freeze it 
will build up at that control volume until the aerody- 
namic forces become larger than the surface tension 
force, which will cause flow into the next control vol- 
ume as runback water. 

2. Runback water flows in the direction opposed 
to the aerodynamic stagnation point. Because of this, 
no runback water can flow into the control volume 
which contains the stagnation point. 

3. Any runback water which the stagnation point 
control volume generates is equally divided into the 
two control volumes on each side of this point. 

Because of the above assumptions, the mass 
balance must be applied to the stagnation point con- 
trol volume first, as the amount of runback water 
present elsewhere is dependent upon the amount 
generated at the stagnation point. The general form 
of the mass balance in terms of the mass fluxes is: 

+ m, b, in = ™e + + >»rb. on, + ™sh + < 38 ) 


The mass flux of impinging water is given by: 

rh K = p (LWC)V„ (39) 

The collection efficiency (P) given above is sim- 
ply the fraction of water droplets in a volume of water 
which impinges upon the surface, meaning that the 
droplet strikes the surface and remains there and 
does not splatter. The collection efficiency can be 
obtained by correlations with the ambient conditions, 
or by a separate computer program which calculates 
the trajectory of the droplets in the air and deter- 
mines if they will strike the surface and if so, if they 
will splatter or remain upon the surface. The collec- 
tion efficiency is calculated by trajectory module in 
LEWICE 2.2. 

The mass lost by water shedding is determined 
independently from the other mass balance terms 
and is based on the Weber number. Water shedding 
can occur from runback watter being sheared from a 
horn or by removal of existing surface water by 
shashing of incoming drops. The amount of water 
shedding is assumed to encompass both water 
shedding and water splashing as both are consid- 
ered to be controlled by the Weber number. The 
empirical expression for the amount lost is given by 

. ( W e~ W e A 

m sh = w >,( — jjjr “ j < 4 °) 

where W e is the Weber number of the runback 
bead and W e c is the critical Weber number to initiate 
shedding. A value of W e c = 500 is currently used. If 
the Weber number is less than W e c then no mass is 
lost by water shedding. 

The term m r is the amount of unfrozen water 

which is not allowed to leave the control volume due 
to surface tension (Weber number) effects. The run- 
back velocity is assumed to be 1/2 of the local air 
velocity at the edge of the boundary layer. The value 
of 1/2 was determined by assuming the water flow is 
shear driven. If the mass flux of water computed to 
runback is more than the flux computed using this 
velocity, the remaining water is kept at that control 
volume until the next time step. 

The mass flux of evaporation is given by: 
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One of the purposes of equations in the previous 
section was to obtain the freezing fraction, N f . This is 
the fraction of impinging water which freezes. Hence, 
the amount of runback water leaving the control vol- 
ume is obtained by multiplying the total amount of 
incoming water by the fraction which does not freeze. 
Hence: 

+ ™rb, ou, = «■ + ™rb. ,•«)( 1 - Nf) (42) 

The mass flux of runback water coming in, 
m rb ( „, is zero at the stagnation point and is calcu- 
lated from the previous volume’s mass balance else- 
where. This is why it is necessary to start the 
calculations at the stagnation point. The right-hand 
side of this equation is compared to the maximum 
runback mass flux, piVg/2. If the amount of available 
water is more than the maximum, the excess water is 
stored in m r until the next time step. The mass flux of 

ice accretion is obtained by solving Eq. 38 for m a and 
substituting the above expressions to obtain: 

m a = (m w - m rb jn )N f - - rii sh (43) 

Ice Shedding 

The assumptions made in the development of an 
ice shedding model are: 

1 . The flow is two-dimensional, that is, the aero- 
dynamic force has x and y-components only. 

2. The forces holding the ice to the surface are 
determined from the bonding strength of the ice. 

3. The bonding strength is dependant only upon 
the temperature at the ice/airfoil interface. 

4. The ice can shed as a whole or in sections. 

5. Each spanwise section acts independently of 
the other sections. 

6. The ice will shed when the net average exter- 
nal forces exceed the net average force holding the 
ice to the surface. The aerodynamic force is com- 
puted from the flow solution as either a macroscopic 
force over the whole airfoil or at each individual loca- 
tion, depending on the shedding mechanism chosen 
by the user. 

The equations relating bonding strength of the 
ice versus temperature were obtained by curve-fitting 

experimental data provided by Scavuzzo 8 9 . The com- 
putational domain takes the x-direction to be along 


the surface of the airfoil in the chordwise direction. 
The y-direction is normal to the airfoil, and the z- 
direction is along the surface of the airfoil in the span- 
wise direction. Shedding ice as a whole assumes 
that the tensile strength of ice is very high. This 
means that if part of the ice adheres to the surface, 
all of the ice will remain. Shedding by node assumes 
that ice has no tensile strength. This means that if the 
adhesive strength at that location goes to zero (the 
ice melts) then it can not be held in place by the 
neighboring nodes. Prediction of surface tempera- 
tures and shedding times from the experimental data 
set indicate that the second assumption is more likely 
for the conditions tested. 

VI. Bleed Air Results 

The bleed air design used for this analysis is 
shown in Fig. 2. It was designed by Jim Van Fossen 
and Phil Poinsette from the NASA Glenn Heat Trans- 
fer Branch based upon their experience with turbine 
blade cooling. Some examples of their previous work 
have been provided for reference 10 ' 13 . The model 
was designed for use as a validation model and was 
not intended to represent any existing design used by 
industry. 

The exterior surface is a NACA0012 airfoil, 
shown in Fig. 3. The interior surface has a uniform 
channel with small holes at the chordwise locations 
shown in Fig. 2. The holes were 0.045” in diameter 
and were spaced 10 diameters, or 0.45” apart. A 
benchtest was performed using liquid crystals to 
determine the heat transfer coefficients on the inside 
surface 14 . A sample heat transfer curve from this test 
is shown in Fig. 4. This model will be tested in FY02 
in the IRT. An example output for an evaporative con- 
dition from the test matrix is shown in Fig. 5. 

A second example case uses five 2D cross sec- 
tions of an engine inlet as the body geometries. The 
first of these cross-sections is shown in Fig. 6. This 

model was taken from a LEWICE3D 15 example case 
provided by Colin Bidwell.The surface pressure coef- 
ficients, collection efficiencies, and external heat 
transfer coefficients were supplied from a LEWICE3D 
output file at each 2D cross-section shown. The inter- 
nal heat transfer coefficients for this case were 
selected such that residual ice would form. A plot of 
these values are shown in Fig. 7. These values were 
then used for each cross-section. Since the inputs 
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and results for each cross-section are very similar, 
only the inputs and results for the first section are 
shown in this paper. The LEWICE 2.2 draft user man- 
ual provides additional analysis of this case. Figure 8 
shows surface temperature predictions for the first 
section. The temperatures and the resultant residual 
ice, shown in Fig. 9, can then be loaded back into 
LEWICE3D for the generation of a simulated 3D 
residual ice shape. The simulated 3D ice shape is 
shown in Figure 10. The resolution of the ice shape in 
this figure can be improved by the addition of more 
streamlines to the model or by interpolation from the 
existing results. Since LEWICE 2.2 also contains a 
prediction of rivulet formation in its model, it should 
be possible in the future to display that information as 
well. 

VII. Electrothermal Validation 

The test article used to generate the electrother- 
mal deicer database was a NACA0012 airfoil with a 
72 inch span and a 36 inch chord. The leading edge 
had seven independently controllable heater zones 
as shown in Fig. 1 1 . The material properties of the 
heater mat are given in Table 1 . Heater zone A (part- 
ing strip) was on continuously for all the runs. The 
other heaters were cycled using various power set- 
tings and on/off times. Results from this test were 
previously reported 15 . 

Development of the test matrix involved the 
selection of two different types of parameters: icing 
parameters (T 0 , LWC, MVD, etc.) and electrothermal 
ice protection system parameters (heater power level 
and heater zone on/off time). The combination of 
both sets of parameters resulted in an extremely 
large number of possible test parameter combina- 
tions. Consequently, the test was restricted to a few 
variations in icing conditions and a wider variety of 
ice protection parameters. Including repeat runs, 
there were 1 1 3 runs in this test entry. 

For each data point, three LEWICE 2.2 outputs 
were generated. In the first set of cases, the program 
used the standard external heat transfer coefficient 
which is calculated for icing runs. A second set of 
outputs were generated using a laminar heat transfer 
coefficient on the external boundary. This assumption 
argues that the airfoil surface is kept relatively clean 
by the de-icer, therefore the boundary layer will 
remain laminar within the impingement limits (natural 


transition normally occurs further downstream). The 
third set of comparisons used the laminar boundary 
layer assumption as well as an assumption of zero 
tensile strength for the ice, allowing each node of ice 
to shed independently of the other. The first two sets 
of output assumed that the ice shed in sections. All of 
these cases have been compared to the LEWICE 2.2 
model. 

Figure 12 shows the comparison on the parting 
strip for Run17 from the test matrix using the stan- 
dard boundary layer transition model. This plot shows 
the comparison for an RTD (resistance temperature 
device) located below the heater. Figures 13 and 14 
show results comparing thermocouples at the sur- 
face and on the inside surface of the airfoil. Figures 
15-20 show the comparisons for sections “B” and “D” 
for the same run. The quality of this comparison is 
typical for the deicing cases in this matrix. Since the 
model is symmetric, only data from one side needs to 
be shown. The fluctuations in the experimental 
results for section "A” occur when section “B” and “C” 
heaters are on. This effect does not occur numeri- 
cally due to the underprediction of those tempera- 
tures. The conditions for this case, including heater 
power and cycle times, is shown in Table 2. 

This comparison shows that section “A” is well 
predicted, while the other zones are underpredicted. 
This result occurred because section “A” is in the 
laminar region of the boundary layer while the other 
zones see the effects of the higher turbulent heat 
transfer coefficient. Therefore, the case was run 
using a laminar boundary layer assumption. The 
results for this case are shown in Figs. 21-29. Fig- 
ures 21-23 show the output for section “A”. Although 
the boundary layer is laminar for the previous case as 
well, this case has a higher temperature as less heat 
is lost to the other sections. Figures 24-29 show a 
larger increase in temperature due to the change in 
assumption. Although the temperatures are slightly 
overpredicted, the comparison has improved signifi- 
cantly. Once again, this result is typical for a large 
majority of de-icing cases. A minority of cases still 
show an underprediction even with the laminar 
assumption. 

The last set of comparisons uses an alternative 
shedding mode than the previous runs. The previous 
runs allowed ice to shed only section by section, 
meaning that any ice which would accumulate could 
only shed if the aerodynamic force exceeded the 
adhesion force over the entire section. The sectional 
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shedding model is a new feature not modeled by pre- 
vious deicer programs. LEWICE/Thermal could shed 
ice an a whole or by each node acting independently. 
This latter assumption was used in the third set of 
comparisons in order to assess its value. For brevity, 
only the heater temperatures are shown. Figures 30- 
32 show the results of using this assumption. While 
the comparisons are similar to those obtained with 
the sectional shedding model, the latter heater cycles 
show a tendency of the nodal shedding cases to 
increase the peak temperature. Once again, these 
results are fairly typical for cases ran using this 
assumption. 

VIII. Quantitative Analysis 

The previous section used a qualitative analysis 
to compare outputs from LEWICE 2.2 to the experi- 
mental data. This section will describe a methodol- 
ogy which was created to quantify those findings. 
The first step in this process was to create a version 
of the software that runs in batch mode. In this mode, 
interactive inputs were automatically answered so 
that the user did not need to be at the display enter- 
ing filenames. This mode is not recommended for 
standard use, as the warning and error messages 
displayed often point to input file errors. Second, a 
script was written to run each of the 113 cases 
sequentially, along with a small data reduction pro- 
gram that converted temperatures to Fahrenheit and 
reduced the output to one value per second to match 
the experimental data. 

The third step involved writing a Visual Basic® 
program in Microsoft Excel®. In this step, each sec- 
ond of experimental data was compared to each sec- 
ond of computational output, with the absolute 
difference and absolute percent differences calcu- 
lated for each comparison. The macro also produced 
the comparison plots shown earlier. These steps 
greatly reduced the man-hours necessary to produce 
a comprehensive comparison which allowed the 
exploration of alternative assumptions shown previ- 
ously. The average differences and percentage differ- 
ences for each thermocouple reading were then 
averaged for that case and finally the results for each 
case were averaged. The differences and percent dif- 
ferences were calculated by the following equations: 

abs difference = |T exp -T LEmcE \ (44) 


l . a-et 100|r p - T LEW/Ce\ 

abs percent difference = — — J 1 — (45) 

1 ex P max ~ ‘ e *P mb , 

Table 3 shows a summary of results from the 
quantitative analysis. There are 69 cases which 
reside under the de-icing category. While the overall 
results are similar for the three sets of comparisons, 
it should be noted that both laminar comparisons 
resulted in an overprediction of temperatures while 
the turbulent cases resulted in an underprediction of 
temperatures. This result indicates that it would be 
possible to produce a calibrated external heat trans- 
fer coefficient which would more accurately predict 
the experimental temperatures. This calibration was 
not performed, as the purpose of the validation was 
to show the comparison capabilities using current 
capabilities only. This table also shows that the per- 
centage difference between LEWICE 2.2 and experi- 
ment is not a function of the ambient temperature 
and that a larger temperature difference is obtained 
for colder conditions which have more power applied 
to them. This result is consistent with the conclusion 
that a calibrated heat transfer coefficient could be 
derived from the experimental results. 

There were 44 cases which had power supplied 
continuously to all of the heaters. Of these cases, 1 7 
were run with a high enough power setting to evapo- 
rate all of the incoming water. The other 27 cases 
had runback ice past the heater zones. Table 3 also 
shows the predictive capability of LEWICE 2.2 for 
these cases. Although the average absolute differ- 
ence for these cases was extremely high, the aver- 
age absolute percentage difference using the laminar 
assumption is similar to the de-icing cases. In this 
case however, a calibrated heat transfer coefficient 
can not be derived. The laminar assumption with sec- 
tional shedding cases resulted in exactly 50% of the 
44 cases being underpredicted and 50% over pre- 
dicted. The nodal shedding assumption increased 
the number over predicted to 53%. For the turbulent 
assumption, 23% were over predicted with 77% 
underpredicted. There may be several factors which 
could explain both the poor prediction of the temper- 
atures as well as their inconsistency. 

The main cause appears to be that while the 
LEWICE 2.2 results were ran with a constant power 
setting, none of the experimental cases were ran in 
that mode. The power was continually adjusted both 
manually and with a temperature controller to obtain 
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a desired temperature response. Figure 33 shows 
the power input for Section A of Run 35A as a typical 
example. Experimentally, the power varied widely for 
the first part of the run, then was stabilized except for 
a period of time where it appears the power was acci- 
dentally turned off. LEWICE 2.2 cannot model this 
type of heater input. All of the LEWICE 2.2 cases 
were ran with a constant power setting at the level at 
which the power stabilized. The comparison was then 
made between experiment and LEWICE 2.2 only in 
the stabilized region. Despite these efforts, the tem- 
peratures were not well predicted for most of the 
evaporative cases. 

IX. Ice Shape Validation 

A quantitative analysis similar to the LEWICE 2.0 
validation effort was also performed with this version 
of the software. A script similar to the one used for 
de-icing comparison was generated to reduce the 
time necessary for this effort. In this effort, the terms 
calculated are the same as those described in the 

validation report 4 with two exceptions. First, the 
upper and lower horn angles are now calculated rela- 
tive to the leading edge of the clean airfoil and not to 
the center of the inscribed circular cylinder. The new 
definition is illustrated in Fig. 35. Second, the ice 
area is calculated by the following equation: 

Area = (fy dx) -(fv dx ) (46) 

J ' iced J ' clean 

In order to compare these results to the experi- 
mental database and to LEWICE 2.0, these parame- 
ters were determined from the previous results. This 
operation was performed using similar scripts to 
those generated earlier. The results of this compari- 
son are shown in Table 4. While the overall compari- 
son of LEWICE 2.0 and LEWICE 2.2 results to 
experiment produce virtually identical results as 
shown in this table, individual cases can produce 
somewhat different ice shapes. Additionally, it is 
worthwhile mentioning that the automated process 
used here produced different values for each param- 
eter when it was ran using LEWICE 2.0 than the val- 
ues obtained from the LEWICE 2.0 validation report. 
Therefore, the automated process may not be an 
accurate assessment of the software’s capabilities for 
ice shape generation. 


X. Conclusions 

Results to date have shown a reasonable com- 
parison with deicing cases which have been ana- 
lyzed. The results show that it would be possible to 
create a calibrated heat transfer coefficient which 
would better predict the experimental data. Evapora- 
tive cases are not well predicted with either the lami- 
nar or turbulent assumption. The turbulent cases 
have shown significant underprediction of tempera- 
ture while the laminar cases show an overprediction 
in half of the cases. Ice shapes generated with 
LEWICE 2.2 have been shown to be quantitatively 
similar to the LEWICE 2.0 results, although the find- 
ings to date have been limited in scope. 
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FIGURE 1. Representation of a Simple Deicer Model 
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FIGURE 3. NACA0012 Airfoil 
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FIGURE 6. Engine Inlet Cross-Section 

FIGURE 4. Experimental Heat Transfer Coefficients for 
Bleed Air Model 



FIGURE 5. Predicted Temperature Distribution for 
Evaporative Bleed Air System 



FIGURE 8. Surface Temperature Distribution for Engine 
Inlet Cross-section 
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FIGURE 9. Residual Ice Shape for Engine Inlet Cross- 
section 



FIGURE 10. Simulated 3D Ice Accretion for Engine Inlet 
Half-Model 


Top View 
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FIGURE 11 . Schematic of Heater Zones 
TABLE 1. Properties and Thicknesses of Deicer Pad 


# 

Material 

Thickness 

(mm) 

Conductivity 

(W/m/K) 

Diffusivity 

(m 2 /s*10 7 ) 

1 

Foam 

3.43 

0.121 

1.65 

2 

Fiberglass 

0.89 

0.294 

1.04 

3 

Elastomer 

0.28 

0.256 

1.50 

4 

Heater 

0.013 

41.0 

120.0 

5 

Elastomer 

0.28 

0.256 

1.50 

6 

Ab. shield 

0.20 

16.3 

40.6 


Run 17 Top, Soctlon A: Haator Tamparatura 
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FIGURE 12. Turbulent Parting Strip Heater Temperature 
Comparison for Run17 
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Run 17 Top, Section A: Outor Surface Temperature 


Run 17 Top, Section B: Outer Surface Temperature 
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FIGURE 13. Turbulent Parting Strip Surface Temperature 
Comparison for Run17 


FIGURE 16. Turbulent Section “B” Surface Temperature 
Comparison for Run17 


Run 17 Top. Section A: Inner Surface Temperature 



Time (a) 

# Section A Inner Surface Temperature (Experimental) [*F] -Section A Inner Surface Temperature ILEWICE) [*F] 


Run 17 Top, Section B: Inner Surface Temperature 



FIGURE 14. Turbulent Parting Strip Substrate 

Temperature Comparison for Run17 


Run 17 Top, Section B: Heater Temperature 



Time (s) 

• Section B Heater Temperature (Experimental) [*F) Section B Healer Temperature (LEWICE) f*F] 


FIGURE 17. Turbulent Section “B” Substrate Temperature 
Comparison for Run 17 


Run 17 Top, Section D: Heater Temperature 



Time (e) 


• Section 0 Heater Temperature (Experimental) [*F] Section D Heater Temperature (LEWICE) [*F] 


FIGURE 15. Turbulent Section “B” Heater Temperature 
Comparison for Run17 


FIGURE 18. Turbulent Section “D” Heater Temperature 
Comparison for Run 17 
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Run 17 Top. Section D: Outor Surface Temperature 
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FIGURE 1 9. Turbulent Section “D” Surface Temperature 
Comparison for Run17 


Run 17 Top. Section D: Inner Surface Temperature 



FIGURE 20. Turbulent Section “D” Substrate Temperature 
Comparison for Run17 

TABLE 2. Conditions for Electrothermal Run 17 
MVD= 20. pm 
c = 0.9144 m 
a = -0.5° 

V°° = 44.7 m/s 
LWC= 0.78 g/m3 
Too = 265.50 K 
Poo = 1 00 kPa 
rh =100% 

Htr. A 7.9 kW/m2 ON always 

Htr. B 13.4 kW/m2 ON 10 sec. OFF 50 sec. 

Htr. C 13.9 kW/m2 ON 10 sec. OFF 50 sec. 

Htr. D 13.5 kW/m2 ON 10 sec. OFF 50 sec. 

Htr. E 13.5 kW/m2 ON 10 sec. OFF 50 sec. 

Htr. F 13.9 kW/m2 ON 10 sec. OFF 50 sec. 

Htr. G 13.6 kW/m2 ON 10 sec. OFF 50 sec. 

Htrs. D-G turned ON after Htrs. B&C turned OFF 
in the cycle. 


Run 17 Top. Section A: Hoater Temperature 
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FIGURE 21. Laminar Parting Strip Heater Temperature 
Comparison for Run17 


Run 17 Top, Section A: Outer Surface Temperature 



Time (a) 

# Section A. Outer Surface Temperature (Expertmental) [*F) Section A Outer Surface Temperature (LEWICE) (*F| 


FIGURE 22. Laminar Parting Strip Surface Temperature 
Comparison for Run17 


Run 17 Top. Soctlon A: Innor Surface Temperature 



FIGURE 23. Laminar Parting Strip Substrate Temperature 
Comparison for Run17 
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Run 17 Top. Soctlon B: Haatar Tomporaturo 


Run 17 Top. Soctlon 0: Hoator Tomporaturo 



Tima (a) 

l Section B Haalar Temperature (Expenmentall |*F] ■■ Section B Healer Temperature (LEWICE) |*F] 



Tima (a) 

• Sactton 0 Heater Tamparatura (Experimental) [*F] Soctlon 0 Haatar Tamparatura (LEWICE) t*F] 


FIGURE 24. Laminar Section “B” Heater Temperature 
Comparison for Run17 


FIGURE 27. Laminar Section “D” Heater Temperature 
Comparison for Run17 


Run 17 Top. Soctlon B: Outor Surface Tomporaturo 


Run 17 Top. Soctlon D: Outor Surface Tomporaturo 
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FIGURE 25. Laminar Section “B” Surface Temperature 
Comparison for Run17 


FIGURE 28. Laminar Section “D” Surface Temperature 
Comparison for Run17 


Run 17 Top. Soctlon B: Innor Surface Tomporaturo 
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FIGURE 26. Laminar Section “B” Substrate Temperature 
Comparison for Run17 


Run 17 Top, Soctlon D: Innor Surfaco Tomporaturo 
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FIGURE 29. Laminar Section “D” Substrate Temperature 
Comparison for Run 17 
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Run 17 Top. Soclion A: Hooter Temperature 



Time (•) 

• Sod loo A Hosier Temperature (Experimental) |*F] Section A Healer Temperature HEWICE I |*F| 


FIGURE 30. Laminar Section “A” Heater Temperature 
Comparison with Nodal Shedding for Run17 


Run 17 Top. Section B: Heater Temperature 
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FIGURE 31. Laminar Section “B” Heater Temperature 
Comparison with Nodal Shedding for Run17 


Run 17 Top. Section D: Heater Temperature 



FIGURE 32. Laminar Section “D” Heater Temperature 
Comparison with Nodal Shedding for Run17 


TABLE 3. Quantitative Comparison of Experimental and 
Computational Deicing Results 



Laminar, 

Sectional 

Shedding 

Laminar, 

Nodal 

Shedding 

Turbulent, 

Sectional 

Shedding 

Overall 

Deicing 

8.6°F, 

25.7% 

9.2°F, 

29.6% 

8.8°F, 

27.2% 

20°F Cases 

6.5 °F, 
26.3% 

7.3°F, 

29.6% 

5.9°F, 

22.4% 

0°F Cases 

1 1 .9°F, 
24.5% 

12.2°F, 

29.6% 

13.2°F, 

29% 

Evap. Anti- 
Icing 

35 °F, 
22.7% 

55 °F, 
35.6% 

88.6°F, 

56% 

Running 

Wet 

22.2°F, 

38% 

25.6°F, 

44.2% 

33.4°F, 

53.9% 





FIGURE 33. Heater Power for Run 35A, Section A 


Run 35A Top. Section A: Heater Temperature 
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FIGURE 34. Temperature Comparison for Run 35A, 
Section A using Laminar Heat Transfer 
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FIGURE 35. Definition of Horn Angle 


TABLE 4.Comparison Of Icing Parameters for LEWICE 2.0 
and LEWICE 2.2 


Parameter 

LEWICE 2.2 
(Automated 
Process) 

LEWICE 2.0 
(Automated 
Process) 

LEWICE 2.0 

(Validation 

Report) 

[S8] Lower 
Icing Limit 
Difference 
(%c) 

4.8 

4.8 

6.1 

[S7] Upper 
Icing Limit 
Difference 
(%c) 

1.35 

1.27 

1.65 

[S6] Lead- 
ing Edge 
Minimum 
Difference 
(% of total 
thickness) 

7.5 

8.1 

5.7 

[S5] Iced 
Area Differ- 
ence (% of 
total thick- 
ness * 
width) 

10.5 

10 

9.9* 

[S4] Lower 
Horn Angle 
Difference 
(degrees) 

30.7 

29.7 

29.6* 


TABLE 4.Comparison Of Icing Parameters for LEWICE 2.0 
and LEWICE 2.2 


Parameter 

LEWICE 2.2 
(Automated 
Process) 

LEWICE 2.0 
(Automated 
Process) 

LEWICE 2.0 

(Validation 

Report) 

[S3] Lower 
Horn Thick- 
ness Differ- 
ence (% of 
total thick- 
ness) 

14.1 

14 

11.7 

[S2] Upper 
Horn Angle 
Difference 
(degrees) 

22.3 

20.9 

16.4* 

[SI] Upper 
Horn Thick- 
ness Differ- 
ence (% of 
total thick- 
ness) 

14.7 

13 

9.7 


* The LEWICE 2.0 validation report used differ- 
ent equations to calculate these parameters, there- 
fore the values cannot be directly compared. 
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